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Abstract 


Previous  studies  have  shown  that  the  Space  Shuttle 
Orbiter  can  achieve  larger  orbit  inclination  changes  using 
an  aerodynamic  turn  than  can  be  obtained  using  a  rocket 
motor  burn.  This  analysis  determines  the  angle  of  attack  and 
bank  angle  histories  which  maximize  the  change  in  inclin¬ 
ation  angle  while  satisfying  final  altitude  and  velocity 
loss  constraints.  The  angle  of  attack  and  bank  angle  are 
modelled  as  polynomial  functions  of  time  with  unknown 
coefficients.  The  optimum  values  of  the  coefficients  are 
determined  by  a  gradient  optimization  technique.  Additionally, 
the  sensitivity  of  the  change  in  inclination  angle  to  changes 
in  the  orbit  apogee  altitude  is  examined.  It  is  shown  that 
the  Space  Shuttle  Orbiter  can  obtain  higher  inclination 
angle  changes  from  orbits  with  higher  apogee  altitudes. 
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USE  OF  AN  AERODYNAMIC  TURN  TO 


MAXIMIZE  THE  ORBIT  INCLINATION 
CHANGE  FOR  THE  SPACE  SHUTTLE  ORBITER 

I .  Introduction 

Background 

When  the  Space  Shuttle  Orbiter  becomes  operational , 
new  spacecraft  capabilities  will  become  available  to  space 
mission  planners.  Not  only  is  the  orbiter  reusable,  but 
for  the  first  time  mission  planners  will  also  have  a 
spacecraft  with  significant  aerodynamic  capabilities. 

In  certain  regimes,  the  orbiter  has  a  lift  to  drag  ratio 
near  two.  When  operating  in  this  regime,  the  shuttle  has 
the  option  of  orienting  the  lift  vector  so  that  a  change 
in  orbit  inclination  is  now  possible  via  an  aerodynamic 

turn  rather  than  the  conventional  rocket  motor  burn  method. 

Wnat  this  offers  the  mission  planner  is  a  twofold  capability. 

One,  the  orbiter  can  use  the  aerodynamic  turn  to  increase 
the  number  of  maneuvers  possible  during  the  course  of  a 
mission  which  would  allow  the  servicing  of  more  satellites. 

Secondly,  the  shuttle  may  trade  fuel  for  increased  payload 
at  launch  and  still  be  able  to  accomplish  a  required 
minimum  number  of  orbital  maneuvers  using  the  aerodynamic 
turn. 

It  is  desirable  to  use  the  aerodynamic  turn  to  achieve 
an  inclination  change  if  the  maneuver  is  more  efficient 
than  the  conventional  rocket  motor  technique.  The 
aerodynamic  turn  is  more  efficient  than  the  rocket  burn 
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if  one  of  the  following  two  conditions  are  met:  (1)  either 
the  aerodynamic  turn  requires  less  velocity  loss  due  to 
air  drag  (Av)  to  achieve  the  same  inclination  change  or 
(2)  more  of  an  inclination  change  is  available  for  the 
same  velocity  loss  when  comparing  the  aerodynamic  turn 
to  the  rocket  burn. 

If  the  aerodynamic  turn  is  more  efficient  than  the 
rocket  burn,  the  orbiter  will  have  several  new  options. 

First,  the  shuttle  can  use  the  aerodynamic  turn  to  achieve 
an  inclination  change  and  then  use  the  rocket  motor  to 
regain  the  velocity  lost  to  aerodynamic  drag.  The  end 
result  would  be  a  larger  orbit  inclination  change  than  would 
be  possible  using  only  the  rocket  motor.  Using  this  method, 
the  'shuttle  orbit  would  have  the  same  shape,  but  a  new 
Inclination.  Secondly,  the  shuttle  could  use  both  the 
aerodynamic  turn  and  the  rocket  motor  to  achieve  a  larger 
inclination  change  than  is  possible  by  either  method  alone. 

The  resulting  orbit  would  have  not  only  a  new  inclination, 
but  a  different  shape.  Obviously,  many  variations  of  these 
two  extremes  are  also  available. 

Previous  research  (Ref  5)  shows  that  for  perigee 
altitudes  of  85  km  and  lower,  the  aerodynamic  turn  maneuver 
does  achieve  the  same  inclination  change  for  less  Av 
cost  than  the  rocket  motor  burn.  This  earlier  study  on 
the  aerodynamic  turn  attempted  to  optimize  the  angle  of 
attack  and  bank  angle  control  history  which  minimized 
the  work  done  by  air  drag  (Av)  while  satisfying  end  conditions 
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that  specified  final  altitude  and  orbit  inclination  change. 
Because  of  this  approach  to  the  problem,  the  optimization 
technique  did  not  function  properly  and  a  non-optimum 
approach  to  the  efficiency  of  the  aerodynamic  turn  was 
finally  taken.  The  aerodynamic  turn  maneuver  was  found 
to  be  more  efficient  than  the  rocket  motor  burn  and  the 
maximum  orbit  inclination  change  obtained  by  the  aerodynamic 
turn  using  a  constant  angle  of  attack  and  bank  angle 
throughout  the  maneuver  was  .78°.  This  earlier  study 
recommended  that  research  be  continued  on  the  subject 
with  a  restructured  optimization  problem  and  that  the 
aerodynamic  turn  maneuver  be  limited  to  95  km  and  lower 
perigee  altitudes. 

Problem  Statement  and  Scope 

This  current  investigation  seeks,  first,  to  find 
another  approach  to  the  optimization  problem  so  that  the 
difficulties  encountered  in  the  previous  research  (Ref  5: 
49-50)  are  eliminated.  If  the  problem  can  be  successfully 
restated  so  that  optimization  is  possible,  the  maximum 
orbit  inclination  change  obtainable  using  the  aerodynamic 
turn  maneuver  is  desired. 

The  aerodynamic  turn  optimization  problem  is  therefore 
restated  as  follows:  Find  the  angle  of  attack  and  bank 
angle  histories  which  maximize  inclination  change  while 
satisfying  end  conditions  that  require  returning  the 
shuttle  to  a  given  final  altitude  with  a  specified  accept¬ 
able  velocity  loss  due  to  air  drag. 
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Because  the  shuttle  must  obey  the  appropriate  equations 
of  motion  during  the  aerodynamic  turn  maneuver,  the  max¬ 
imization  of  inclination  change  is  an  optimal  control  problem. 
While  there  are  several  ways  to  solve  an  optimal  control 
problem,  the  one  chosen  here  is  the  Second-Order  Control 
Parameter  technique  (Ref  7).  In  this  method  the  orbiter 
controls  are  represented  suboptimally  as  polynomial 
functions  of  time.  The  coefficients  of  the  polynomials 
are  chosen  to  optimize  the  performance  index  (maximum 
inclination  change)  and  satisfy  the  end  conditions  (final 
altitude  and  velocity).  Higher  order  control  polynomials 
are  investigated  to  determine  their  effect  on  inclination 
change  and  what  order  controls  are  necessary  to  accurately 
approximate  the  optimal  controls.  Additionally,  a  sensit¬ 
ivity  analysis  is  conducted  on  various  apogee  altitudes 
to  determine  their  effect  on  inclination  change. 

The  remainder  of  this  paper  is  arranged  in  the  following 
manner:  Chapter  II  discusses  the  orbiter  equations  of 

motion  and  the  model  atmosphere  used  in  this  study. 
Explanations  of  the  methods  used  to  generate  coefficients 
of  lift  and  drag  are  also  given.  Chapter  III  covers  the 
optimization  routine  used,  why  it  was  selected  and  how  it 
is  implemented.  Chapter  IV  reports  on  how  the  initial 
values  of  control  coefficients  and  end  conditions  used  in 
the  optimization  routine  are  selected.  Chapter  V  lists 
the  results  and  conclusions  generated  by  this  study. 


II. 


System  Dynamics  and  Model  Atmosphere 


Background 

The  Space  Shuttle  Orbiter  Is  a  dynamic  system  and 
must,  therefore,  obey  Newton’s  laws  of  motion.  Specifically, 
Newton's  second  law  relates  the  aerodynamic  and  gravitational 
forces  acting  on  the  orbiter  to  the  orbiter  acceleration, 
EF=ma.  This  relationship  can  be  expressed  as  a  set  of 
differential  equations  which  describe  the  motion  of  the 
shuttle  as  it  passes  through  the  atmosphere.  B>'  ^ore  these 
equations  of  motion  can  be  derived,  a  coordinate  reference 
frame  must  be  specified. 

Coordinate  Reference  Frame 

For  this  study,  two  coordinate  reference  frames 
are  specified.  The  first,  an  X,  Y,  Z  frame  is  earth  centered, 
non-rotating,  and  inertial  (Fig  1).  The  X,  Y  axes  are 
in  the  earth  equitorial  plane  and  Z  is  chosen  to  complete 
a  right-hand  coordinate  frame. 

The  second  frame  is  the  V  M,  L  which  is  fixed 

at  the  center  of  mass  of  the  orbiter  (Fig  1).  The 

Vrel  axis  P°ints  into  the  relative  wind  and  the  angle 

between  the  V  ,  axis  and  the  shuttle's  longitudinal  axis 
rel 

is  the  orbiter  angle  of  attack  a.  The  M  axis  establishes 
a  local  horizontal  for  the  shuttle  and  L  is  chosen  to 
complete  the  right-hand  coordinate  frame.  The  angle 
between  the  shuttle  lift  vector  and  the  L  axis  establishes 
the  angle  of  bank  <f> . 
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The  Orbiter's  aerodynamic  forces  are  readily 


expressed  in  terms  of  the  body  fixed  frame.  These  forces 
are 

LIFT  =  i  LIFT  j  ( cos<}> L  -  sin^>M) 

-V  -* 

DRAG  =  |  DRAG  1  V  , 

1  '  rel 

where  | LIFT]  and  | DRAG j  represent  the  scalar  magnitude 

of  the  lift  and  drag  vector,  respectively.  V  , ,  M,  and  L 

->  -»■  rel 

are  unit  vectors  of  the  V  , ,  M,  L  frame.  The  position 

rel5  + 

and  velocity  of  the  orbiter  expressed  in  terms  of  X,  Y,  Z 
frame  unit  vectors  is 


r  =  XI  +  YJ  +  ZK 

v  =  XI  +  YJ  +  ZK 

•  •  • 

where  X,  Y,  Z,  X,  Y,  and  Z  are  the  system  state  variables. 
The  vrelJ  L  vectors  canbe  expressed  as 


V  ,  =  v  +  wxr 
rel 


M 


=  V  ,xr. 
rel 


=  MxV 


rel 

-r 

where  w  =  earth  angular  velocity.  Substitution  of  the  expressions 

-y  -*■  -y  -y  -y 

for  r  and  v  into  these  equations  allows  M,  and  L  to  be 

expressed  in  terms  of  the  state  variables  and  the  inertial 
frame  unit  vectors.  This  information  allows  the  generation 
of  a  transformation  matrix  T,  which  will  convert  forces 
expressed  in  terms  of  the  body  frame  to  forces  expressed 
in  the  inertial  reference  frame  (Ref  5:16). 
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Equations  of  Motion 

In  addition  to  lift  and  drag,  there  is  also  a  gravity 
force  acting  on  the  orbiter.  The  forces  acting  on  the 
shuttle  expressed  in  the  inertial  frame  are: 

IP  =■  GravityXIZ  +  Llftm  +  DragXYZ 

However,  drag  and  lift  forces  are  expressed  in  the  body 

frame  and  will  have  to  be  converted  to  the  inertial  frame 

using  the  transformation  matrix  T.  Therefore,  the  summation 

of  forces  is  rewritten  as : 

-► 

ZF  =  Gravity  +  T(Li?t„  ML)  +  T(Dragy 

rel  rel 

The  aerodynamic  forces  listed  above  can  be  written 
as  three  second-order  differential  equations  (Ref  5:17). 
These  equations  can  be  further  reduced  to  six  first-order 
equations  which  become  the  shuttle  equations  of  motion. 

These  equations  of  motion  are  listed  below. 

Letting  w^  =  X 


w6  =  Z 


the  equations  of  motion  are 
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w 4  =  -(GM/r3)X  +  (  LIFT  /M)T  -  (  DRAG  /M)T 

^  ^  -*■ 
w c  =  -(GM/rJ)Y  +  (  LIFT  /M)T_  -  (  DRAG  /M)T,. 
b  +  3  +  A 

*6  =  -(GM/r3)Z  +  (  LIFT  /MjTg  -  (  DRAG  /M)Tg 

where  T.  i  =  1  to  6  is  the  appropriate  element  of  the 
transformation  matrix  (Ref  5:18). 

Assumptions 

The  equations  of  motion  listed  above  are  based  on  the 
following  assumptions: 

(1)  The  earth  is  spherical. 

(2)  The  earth  is  inertial. 

(3)  The  crbiter  is  a  point  mass.  All  forces  act 
through  this  point  which  is  the  center  of  gravity. 

(4)  No  aerodynamic  sideslip  occurs  in  the  atmosphere. 

(5)  Atmospheric  winds  rotate  with  the  earth. 

It  is  possible  to  assume  the  earth  is  spherical  because  the 
aerodynamic  turn  is  initiated  from  an  equitorial  orbit. 

Final  state  inclination  is  still  nearly  equitorial  and 
the  non-spherical  effects  of  the  earth  are  neglibible. 

The  earth  is  treated  as  inertial  since  the  time  spent 
In  the  atmosphere  during  the  aerodynamic  turn  is  quite  short 
(approximately  10  min)  and  the  non-inertial  effects  of 
the  earth  would  not  be  a  factor  In  this  length  of  time. 

It  is  reasonable  to  expect  that  a  coordinated  turn  will 
be  flown  during  the  aerodynamic  maneuver  so  that  no  moments 
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will  be  generated  and  the  orbiter  can  be  approximated 
as  a  point  source.  Similarly,  the  maneuver  will  be  flown 
so  that  the  relative  velocity  vector  is  nearly  along  the 
longitudinal  axis  of  the  orbiter  in  order  to  keep  the 
heat  generated  by  hypersonic  speeds  through  the  atmosphere 
on  the  main  heat-protective  surfaces.  Heat  generated 
on  the  side  and  top  of  the  orbiter  will  be  kept  to  a 
minimum  if  sideslip  is  minimized.  Finally,  the  upper 
atmosphere  winds  are  much  too  complex  to  model  exactly 
in  a  study  of  this  size.  However,  to  treat  the  winds  at 
altitude  as  inertial  would  be  too  much  of  a  simplification. 
Accordingly,  the  main  effect  of  the  winds  on  the  shuttle 
orbit  are  modelled  by  assuming  the  winds  to  rotate  with 
the  earth  (Ref  5:12-13). 

Lift  and  Drag  Computations 

The  expressions  used  to  generate  Lift  and  Drag 
have  already  been  listed  in  terms  of  the  reference  frames. 
Restated  in  terms  of  the  state  variables,  the  lift  and 
drag  equations  are 

LIFT  =  isp(X,Y,Z)CL(X,Y,Z,X,Y,Z)V^el(X,Y,Z,X,Y,Z)S 

DRAG  =  Jap(X,Y,Z)CD(X,Y,Z,X,Y,Z)V^el(X,Y,Z,X,Y,Z)S 

The  reference  surface  area  (S)  of  the  orbiter  is  a  known 
constant  value  and  the  relative  velocity  can  readily  be 
obtained  from  integration  of  the  equations  of  motion. 
However,  the  values  for  the  atmospheric  density  (p) 
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and  the  coefficients  of  lift  and  drag  (CT  and  Cn)  will 

J-i  D 

have  to  be  calculated. 

The  orbiter  coefficients  of  lift  and  drag  are  listed 
in  Tables  I  -  IV  (Ref  6).  These  coefficients  are  a 
functions  of  angle  of  attack  (alpha)  and  the  viscous 
parameter  (VBAR).  While  angle  of  attack  is  directly 
obtainable,  because  it  will  be  shown  to  be  a  control 
variable,  the  VBAR  term  is  a  function  of  the  orbiter  state. 


VBAR  =  M  (C  /Re) 2 

OO  x  00  7 

where  M  =  freestream  mach  number 

oo 

^(TjM^)  =  proportionality  factor  for  linear-viscosity 
temperature  relationship 

Re(p ,vrel, 3, Kinetic  Viscosity)  -  Freestream  Reynolds 

Number 

From  the  equations  listed  for  Lift,  Drag,  and  VBAR  it  is 
apparent  that  a  means  of  determining  the  state  of  the 
atmosphere  is  needed. 

Model  Atmosphere 

The  model  chosen  to  represent  the  atmosphere  in  this 
study  is  the  1962  Standard  Atmosphere  (Ref  11).  While 
a  more  current  version  of  this  model  does  exist,  the 
1962  data  is  used  because  the  Space  Shuttle  Orbiter  Lift 
and  Drag  coefficients  are  based  on  the  1962  model. 


The  simplifying  assumption  that  geopotential  altitude  is 
equal  to  geometric  altitude  is  also  used  here  and  results 
in  a  five  to  seven  percent  deviation  from  the  1962 
Standard  Atmosphere  Model  (Ref  5:10).  This  assumption 
allows  the  variation  of  gravitational  acceleration  with 
altitude  to  be  neglected. 

With  this’ model,  the  density,  temperature,  molecular 
weight,  and  viscosity  of  the  atmosphere  at  any  altitude 
out  to  700  km  is  directly  available.  These  atmospheric 
parameters  are  used  in  the  computation  of  VBAR  (Ref  5:22). 
With  VBAR  calculated  and  the  angle  of  attack  known, 
orbiter  lift  and  drag  coefficient  tables  can 
be  entered  and  and  can  be  determined  (Tables  I  -  IV). 

A  bivariate  interpolation  scheme  is  used  to  determine 
the  values  for  CL  and  when  VBAR  and  angle  of  attack 
are  between  the  tabulated  values  listed  in  these  tables. 

For  values  of  VBAR  between  .01  and  .08  the  and 
0^  values  in  Tables  I  -  IV  were  determined  experimentally, 
while  for  VBAR  values  between  .08  and  5.2,  CL  and 
values  were  determined  analytically  via  the  Lockheed 
Bridging  Formula  (Appendix  A) .  This  formula  bridges  the 
transition  between  free -molecular  flow  and  continum  flow 
(Ref  5=23). 
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50.0  .469620  .262590  .123790  .061020  .080150 


Table  ill 


.327700  1. 334140  1.346280  1.348370  1.370520  1.580650 


Summary 

The  orbiter  system  dynamics  are  represented  by  differ¬ 
ential  equations  which  form  the  equations  of  motion.  Inte¬ 
gration  of  the  equations  of  motion  along  the  orbit  allows 
the  determination  of  the  shuttle  state  at  each  integration 
time  step.  With  the  state  of  the  shuttle  known,  the  model 
atmosphere  is  used  to  determine  the  atmospheric  parameters 
required  for  calculation  of  orbiter  lift  and  drag.  Lift 
and  drag  values  are  then  used  in  generating  an  updated 
orbiter  state. 

The  orbiter  state  is  expressed  in  terms  of  eight 
variables;  X,  Y,  Z,  X,  Y,  Z,  4 ,  a .  Since  the  equations  of 
motion  are  six  first-order  differential  equations  and 


there  are  eight  unknowns,  two  of  the  unknown  values 
are  controlled.  For  this  study,  angle  of  attack  (a)  and 
angle  of  bank  (<}>)  are  the  controlled  values.  The  resulting 
solution  of  the  system  state  is  then  checked  and  changes 
are  made  to  the  controlled  values  which  result  in  an  improved 
state.  This  process  is  repeated  through  an  optimization 
routine  which  generates  final  values  for  the  two  variables 
that  result  in  the  "best"  system  state. 


Ill .  Optimization  Technique 

Background 

The  purpose  of  this  thesis  is  to  find  the  angle  of 
attack  and  bank  angle  histories  which  obtain  the  maximum 
orbit  inclination  change  during  an  aerodynamic  turn  maneuver. 
Since  the  initial  state  of  the  orbiter  and  functions  of 
the  final  state  are  prescribedj  and  the  orbiter  is  constrained 
by  the  equations  of  motion,  this  problem  can  be  treated 
as  a  classical  optimal  control  problem  where  the  goal  is  to 
minimize  the  performance  index,  I.  A  description  of  the 
elements  of  the  optimal  control  problem  as  it  applies  to  the 
aerodynamic  turn  maneuver  follows. 

Optimal  Control 

In  the  aerodynamic  turn  maneuver  the  performance  index 

to  be  minimized  is  negative  inclination  change  (-Ai). 

Therefore,  the  problem  is  a  maximization  problem.  To  calculate 

the  inclination  change,  the  angular  momentum  vector  of  the 

final  state  (Hf)  is  compared  to  the  angular  momentum  vector 

of  the  initial  state  (H.).  These  vectors  are  chosen  since 

1 

they  involve  the  complete  orbiter  state  and  are  normal  to 
the  orbit  plane  (Fig  2).  The  inclination  change  is  obtained 
by  determining  the  change  in  direction  between  and  . 
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Figure  2.  Initial  and  Final  Angular 
Momentum  Vectors 
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Since 


H 


r  x  v 


where  r  =  orbiter  position  vector 

v  =  orbiter  velocity  vector 

the  inner  product  of  and  fL  gives  the  expression  for  the 
change  in  inclination: 


41  «  cos-1  (rf  x  vf>  (rf  x  V 


|rf  x  vf| 


rf  x  vf| 


where  |r  x  v|  is  the  scalar  magnitude  of  the  vector  cross 
product . 

The  initial  position  and  velocity  of  the  orbiter  is 
specified  in  terms  of  the  state  variables.  The  shuttle 
proceeds  from  these  initial  conditions  to  a  specified  set 
of  end  conditions  while  maximizing  the  inclination  change. 
The  time  required  to  complete  this  maneuver,  tf,  is 
unspecified.  For  this  problem  the  specified  end  conditions, 
M,  are  altitude  and  velocity  at  tf.  The  final  altitude 
is  to  equal  the  initial  altitude  and  final  velocity  is  to 
equal  initial  velocity  minus  an  acceptable  velocity  loss 
due  to  air  drag  (Av) .  Since  the  orbiter  is  a  dynamic 
system,  it  is  subject  to  orbit  constraints  which  are  the 
equations  of  motion  listed  in  Chapter  II.  Also,  the  shuttle 
is  subject  to  the  control  constraints  that  limit  angle  of 
attack  (20°£a  £  50°)  and  bank  angle  (0°  <_  <t>  <_  180°). 


Suboptlmal  Control 

The  optimal  control  problem  as  outlined  above  is 
reduced  to  a  suboptimal  control  problem  by  approximating 
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the  control  variables  by  a  functional  form  which  involves 
a  number  of  unknown  parameters  (Ref  8).  Suboptimal  control 
is  chosen  in  preference  to  optimal  control  because  the 
algorithm  to  solve  the  optimal  control  problem  requires 
controls  to  be  specified  as  a  function  of  time.  The  time 
span  of  the  aerodynamic  turn  could,  therefore ,  generate 
large  control  vectors.  In  addition,  the  values  to  be  guessed 

in  the  optimal  control  algorithm  do  not  have  a  physical  signif¬ 
icance.  In  the  suboptimal  approach,  there  is  not  only  a  very 

small  number  of  control  coefficients  to  be  guessed,  but  they  also 
have  a  physical  relevence  which  makes  the  initial  guesses 

much  easier  to  generate.  For  this  analysis,  angle  of  attack 
and  bank  angle  are 

m 

a  -  E  b .  C . 

1  =  0 

n 

<t>  =  E  c.C, 

i=0  1  1 

where  b1  and  c^  are  the  unknown  parameters  and  Ch  is  the 
functional  form.  Chebyshev  polynomials  are  chosen  as  the 
functional  form  because  they  are  orthogonal.  This  prevents 
the  matrices  used  by  the  optimization  scheme  from  becoming 
singular  due  to  possible  linear  dependence  of  the  polynomials . 

To  simplify  the  problem,  the  Chebyshev  polynomials  should 
be  defined  over  a  constant  interval.  Through  the  use  of 
Long's  transformation,  the  free  final  time  problem  is 
transformed  into  a  fixed  final  time  problem.  Long's 
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transformation  defines  a  non-dimensional  time 
t  =  t/t 

where  the  range  of  t  is  £o , lj .  The  Chebyshev  polynomial 
may  then  be  defined  over  the  interval  M  .  The  final  time  (tf) 
is  now  another  unknown  parameter  of  the  problem.  Thus,  the 
vector  of  unknown  parameters,  a,  is 


Since  the  problem  is  now  treated  suboptimally ,  one 
objective  is  to  select  the  correct  number  of  controls, 
m  and  n,  that  accurately  approximate  the  optimal  controls. 
If  the  order  of  the  control  polynomials  is  too  small, 
the  results  will  not  accurately  represent  the  optimal 
controls.  On  the  other  hand,  if  too  high  an  order  control 
polynomial  is  selected,  the  small  gain  in  the  performance 
index  will  not  warrant  the  large  Increase  in  computer  time 
required  to  obtain  the  converged  solution.  The  order  of 
the  control  polynomial  representing  the  angle  of  attack 
need  not  necessarily  equal  the  order  of  the  control  poly¬ 
nomial  representing  the  bank  angle. 

Second-Order  Parameter  Optimization 

The  optimal  control  problem  may  now  be  restated  as 
a  parameter  optimization  problem,  where  the  performance 
index,  the  end  condition  constraints,  the  control  variable 
inequality  constraints,  and  the  differential  constraints 
all  depend  on  only  the  unknown  parameters.  The  solution 
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1 


to  this  problem  may  be  found  by  forming  the  augmented 
performance  index,  F, 

F(a,v)  =  1(a)  +  vTM(a) 

where v  is  the  vector  of  constant  Lagrange  multipliers. 

The  conditions  to  be  satisfied  by  an  extremal  point  are 

F^(a,v)  =  0  and  M(a)  =  0. 

3. 

The  values  of  the  parameters  and  the  Lagrange  multipliers 
are  now  guessed.  Unless  the  guess  is  the  extremal  value,  it 
will  be  observed  that 

F^  0  and  M  ^  0 . 

.  a 

T 

In  order  to  drive  F  and  M  to  zero,  these  quantities  are  then 

3. 

linearized  about  the  guessed  values  of  a  and  v,  so  that 

6FT  =  F  6a  +  MT5v  and  6M  =  M  6a 
a  aa  a  a 

where 

a  is  a  p  x  1  vector  of  parameters 

s  fc 

F  is  a  1  x  p  vector  of  1  partial  derivitives 

cl 

F  is  a  p  x  p  matrix  of  2nd  partial  derivitives 

did. 

S  t 

Ma  is  a  m  x  p  matrix  of  1  partial  derivitives 
6 ( )  is  the  variation  of  ()  such  that  6 ( )®( )  w~(  )  ^ 

T 

Since  we  desire  (F  or  M)  =  0,  the  following  relations 

a  new  ° 

are  valid 

6FT  =  -PFT  and  6M  =  -QM 
a  a 


where  Q  and  P  are  weighting  factors  which  control  the  end 
condition  satisfaction  and  optimization,  respectively. 


The  expressions  for  <5a  and  6  v  can  now  be  generated  so  that 


- 1  T  T 

6a  =  _f  PF  +  M  6  v) 
aa  a  a 


6  v =  (M  F'1MT)~1  (-PM  F_1FT  +  QM). 
a  aa  a  v  a  aa  a 

Rather  than  guessing  an  initial  value  for  v}  it  can 
be  computed  using  the  gradient  technique  with  F  set 
equal  to  the  identity  matrix  and  6a  set  equal  to  zero; 
then 


v  =  (M  )_1( (Q/P)M-M  ij) . 

Cl  d  CL  d 

If  it  is  desired  to  use  a  gradient  technique  to  compute 
the  changes  in  a,  those  changes  are  given  by 


5a  =  -P(fJ) 

Cl 

Any  reasonable  value  for  the  control  coefficients  can  be 
guessed  as  long  as  the  guess  conforms  to  the  control 
constraints.  However,  in  the  next  chapter  an  approach 
to  guessing  accurate  values  for  the  coefficients  will  be 
explained . 


IV.  Optimization  Starting  Values 

Previous  research  (Ref  5)  indicated  the  atmospheric 
turn  maneuver  is  profitable  only  at  perigee  altitudes  of 
85  km  and  below.  Table  V  lists  the  orbit  inclination 
changes  which  resulted  from  using  a  60°  angle  of  bank  at 
various  perigee  altitudes. 

Table  V 

Inclination  Changes  Obtained  From  Previous  Research 


85  km 

80  km 

* 

75  km 

70  km 

.  22° 

.4  3° 

.55° 

OO 

O 

*  The  70  km  perigee  orbit  uses  a  bank  angle  of  50°. 

The  bank  angle  of  60°  was  chosen  initially  since  it 
results  In  the  most  inclination  change  without  allowing  the 
maneuver  to  become  a  total  reentry  that  will  not  leave  the 
atmosphere . 

Only  perigee  altitudes  of  80  km  are  considered  in  this 
study.  This  choice  is  made  since  inclination  changes 
achieved  in  an  orbit  with  a  perigee  higher  than  80  km  are 
extremely  small  and  the  viscous  parameter,  VBAR,  was  often 
less  than  0.01  at  perigee  altitudes  of  75  km  and  less. 

The  orbiter  lift  and  drag  coefficient  data  available  is 
limited  to  a  VBAR  greater  than  or  equal  to  0.01. 


25 


Orbit  Selection 


Once  a  perigee  altitude  is  established,  the  selection 
of  an  apogee  altitude  will  allow  the  determination  of  the 
orbital  elements.  Since  the  orbit  used  in  the  previous 
study  is  not  specified,  three  different  orbits  are  used 
here.  A  low  apogee  of  300  km,  a  medium  apogee  of  500  km  and 
a  high  apogee  of  650  km  are  used.  The  high  apogee  is 
selected  to  remain  below  the  hazard  posed  to  manned  flights 
by  the  Van  Allen  radiation  fields.  All  orbits  are  considered 
initially  to  be  direct  equitorial. 

The  three  orbits  used  in  this  study  are  initially  at 
a  zero  degree  inclination  and  have  an  80  km  perigee.  The 
mean  earth  radius  used  for  computations  is  6378.135  km.  The 
Earth  gravitational  parameter  (y)  used  for  this  study  is 

c:  o  p 

3.986012  x  10  km  /sec  .  The  following  two-body  relation¬ 
ships  are  used  to  calculate  the  necessary  orbital  elements 
for  each  of  the  three  reference  orbits  and  assume  no  air 


The  results  of  these  calculations  are  listed  in  Table  VI. 
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Table  VI 

Reference  Orbit  Elements 


*  Velocity  listed  is  at  apogee  (v=l8o°) 


Initial  Position  and  Velocity  Selection 

There  is  enough  information  available  now  to  initialize 
the  optimization  routine  with  the  position  and  velocity  of 
the  orbiter  at  apogee  and  integrate  the  entire  orbit. 
However,  it  is  assumed  that  significant  orbit  plane  change 
will  occur  only  in  atmosphere  of  a  certain  density.  Plane 
change  above  this  altitude  will  be  neglible.  From  a  study 
of  the  1962  Standard  Atmospheric  model  (Ref  11),  an  altitude 
of  100  km  above  the  earth's  surface  and  below  is  assumed  as 
the  level  of  sufficient  density  for  measurable  plane  change. 
The  orbiter  is  inside  the  atmosphere,  as  defined  here,  for 
approximately  10  minutes  per  orbit.  It  is  therefore 
desirable  to  pick  a  starting  position  and  velocity  for 
each  reference  orbit  that  is  just  above  the  100  km  altitude. 
This  will  reduce  the  amount  of  computer  time  needed 
each  iteration  by  neglecting  that  portion  of  the  orbit 
in  which  no  appreciable  plane  change  occurs.  To 
determine  this  starting  position  and  velocity,  the  inte¬ 
gration  routine  is  run  from  apogee  with  angle  of  bank  and 
angle  of  attack  set  to  zero  and  the  time  parameter  set  to 
the  period  of  the  orbit.  Integration  by  this  method 
gives  the  position  and  velocity  vector  for  each  time  step 
and  from  this  a  height  above  the  Earth  is  easily  calculated. 
Integration  of  the  orbit  in  this  manner  results  in  a  perigee 
one  to  two  kilometers  lower  than  that  calculated  using  the 
two-body  problem  assuming  no  air  drag.  These  results  are 
still  valid,  however,  since  integration  using  angles  of 
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attack  and  bank  not  equal  to  zero  results  in  a  shuttle 
perigee  between  79  and  82  kilometers,  depending  on  the 
orbit . 

From  these  integrations  of  the  orbit,  initial  position 
and  velocity  vectors  are  chosen.  These  vectors  are  listed 
in  Table  VII  .  The  position  and  velocity  vectors  become 
the  two  constraints  for  the  optimization  routine.  The 
program  is  designed  to  return  the  shuttle  to  the  starting 
height  on  the  other  side  of  perigee  and  to  the  scalar 
value  of  the  starting  velocity  minus  a  specified,  accept¬ 
able  velocity  loss  due  to  air  drag. 

Control  Coefficient  Selection 

It  is  now  necessary  to  initialize  the  optimization 
routine  with  values  for  the  time,  angle  of  attack,  and 
angle  of  bank  control  parameters.  The  optimization  routine 
is  designed  with  weighting  factors  that  determine  the 
priority  the  routine  places  on  either  achieving  the  end 
constraints  or  optimizing  the  plane  change.  Initially, 
emphasis  is  placed  in  meeting  the  end  constraints  until 
they  are  very  close  to  the  specified  values.  At  that  time, 
the  optimization  weighting  function,  which  has  been  kept 
small,  is  slowly  increased.  In  order  to  minimize  the  computer 
time  necessary  to  meet  the  end  conditions,  the  control 
parameters  which  get  the  final  states  closest  to  the  end 
conditions  are  chosen  to  be  the  initial  guess. 
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able  VII 


ORBIT 


ORBIT 


Initial  Position  and  Velocity  in  Inertial  Frame* 


INITIAL  POSITION 


80/300 

-3103.9951  i 

5724.3016  j 

80/500 

-3238.5651  i 

5680.0536  j 

80/650 

-4606.8106  i 

4636.5561  j 

INITIAL  VELOCITY 


80/300 

-6.8503  i 

-3.8470  J 

80/500 

-6.7306  i 

-4.0819  J 

80/650 

-5.7871  i 

-5.4225  J 

*  All  position  and  velocity  vectors  are  initially 
equitorial  and  have  no  k  components. 
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To  find  the  best  starting  parameters,  over  150  inte¬ 


gration  runs  were  made  for  a  single  orbit  with  various  angles 
of  attack  and  bank  angle  values.  This  allows  a  dete^mina*- ion 
of  how  well  various  combinations  of  the  control  parameters 
meet  the  constraints.  A  listing  of  the  pertir.er.t  results 
from  these  runs  fcr  the  80/300,  80/500,  80/650  orbits 
is  listed  in  Tables  VIII, IX,  and  X  respectively. 

The  first  value  listed  in  Tables  VIII,  IX,  and  X. 
is  the  difference  between  the  starting  and  the  final  altitude 
at  the  end  of  integration.  A  positive  value  means  the 
shuttle  finished  higher  than  the  initial  altitude.  The 
second  value  in  the  tables  is  the  difference  in  initial 
velocity  loss  for  the  maneuver.  In  this  case,  the  allowable 
velocity  loss  is  60  mps .  Positive  velocity  values  mean 
the  orbiter  lost  less  than  60  mps  during  the  orbit.  A 
negative  value  means  the  orbiter  lost  not  only  the  60  mps, 
but  also  the  additional  listed  velocity.  The  third  value 
in  the  tables  is  the  inclination  change  that  results  from 
those  control  parameters.  With  this  data,  the  "best" 
possible  control  parameters  can  b'e  chosen  that  will 
initially  allow  the  optimization  routine  tc  not  only  come 
closest  to  the  desired  end  conditions,  but  also  give  the 
largest  possible  inclination  change. 


Table  IX 


End  Condition  Miss  Errors  80/500  Orbit 


ANGLE 

OF 

BANK 


ANGLE  OF  ATTACK 


10 


15 


20 


25 


30 


o 


ALT  MISS 

30.7 

29.3 

22.5 

VEL  MISS 

-.011 

-.016 

-.018 

Li 

0  . 

0  . 

0  . 

ALT  MISS 

31-9 

29.9 

22.7 

VEL  MISS 

-.010 

-.016 

-.018 

A  i 

.03 

.09 

.  0  9 

ALT  MISS 

31.9 

29.1 

2 1  •  9 

-.010 

-  •  0 1 0 

-.018 

Ai 

.  06 

.  07 

.09 

ALT  MISS 

31.0 

27.1 

19 . 5 

VEL  MISS 

-  .  010 

-.015 

-  .017 

Ai 

.08 

.11 

•  13 

ALT  MISS 

30.2 

26 . 6 

16 . 0 

VEL  MISS 

-.009 

-.019 

-.016 

Ai 

.11 

.  19 

.  17 

km 

kcs 


r-pr 
~  o 


km 


kps 

deg 

km 

KpS 

deg 

km 

kps 

deg 

km 

kps 

deg 


Av  =  60  mps 
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Table  X 


End  Condition  Miss  Errors  8C/650  Orbit 


ANGLE 

OF 

BANK 


ANGLE  OF  ATTACK 


10° 


15° 


20 


20°  25° 


ALT  MISS 

159-7 

152.5 

140.3 

VEL  MISS 

-.164 

-.366 

-.363 

Ai 

0. 

0  . 

0 . 

ALT  MISS 

161.1 

152.6 

141.0 

VEL  MISS 

-.164 

- .  16  6 

-.163 

Ai 

-03 

.04 

•  05 

ALT  MISS 

mam 

153.2 

135-9 

VEL  MISS 

hd 

-.165 

-.16  2 

Ai 

.07 

.09 

.10 

ALT  MISS 

157.6 

149.7 

136.3 

VEL  MISS 

-.163 

-.164 

-.161 

Ai 

.10 

•  13 

.  15 

A.LT  MISS 

156.6 

14  7.6 

134.3 

VEL  MISS 

-.163 

-.163 

-.159 

Ai 

.14 

.  17 

.  20 

km 

kps 

deg 

km 

kps 

deg 

km 

kps 

deg 

km 

kps 

deg 

km 

kps 

deg 


Av  =  60  mps 
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V.  Results  and  Conclusions 


Optimization  Analysis 

A  second-order  technique  to  generate  updated  control 
coefficients  is  preferred  to  a  gradient  technique  because 
of  faster  convergence.  In  the  suboptimal  approach  to  the 
aerodynamic  turn  maneuver  that  was  outlined  in  Chapter  III, 
the  second-order  parameter  method  requires  the  computation 
of  the  F  matrix,  a  second  derivative  matrix,  by  numerical 

ctd 

means.  The  numerical  approach  to  solving  for  F  is  chosen 

d.3. 

over  an  analytical  approach  because  it  requires  considerably 
less  computer  time.  However,  numerical  differentiation 
does  require  perturbations  in  the  quantities  the  derivative 
is  taken  with  respect  to.  For  example,  numerical  computation 
of  F  requires  perturbation  of  the  a's. 

cl  3. 

Early  in  the  control  parameter  optimization  search, 

problems  were  encountered  with  the  second  derivatives.  The 

values  of  the  elements  of  the  F  matrix  were  found  to 

aa 

be  changing  erratically  by  a  power  of  ten.  Since  the 

size  of  the  perturbations  of  the  a’s  affect  the  behavior 

of  the  F  matrix,  a  parametric  study  using  different  size 

perturbations  was  conducted  to  determine  whether  or  not  the 

erratic  behavior  of  the  F  matrix  could  be  corrected. 

aa 

After  evaluating  the  effects  of  various  size  pertur¬ 
bations  on  the  F  matrix,  it  was  determined  that  continued 
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efforts  to  optimize  using  the  second-order  parameter 

method  would  fail  to  produce  a  sequence  of  control 
histories  which  would  ultimately  converge.  A  switch 

to  the  gradient  technique  was  made  at  this  time  to 

see  if  the  optimization  routine  would  move  in  a  consistent 

direction  toward  convergence. 

While  the  gradient  technique  did  continually  move 
toward  an  optimum  control  history  and  improve  both  the 
end  conditions  and  inclination  change,  it  was  very  costly 
in  terms  of  computer  time.  The  results  listed  in  this 
chapter  required  over  32000  CP  seconds  of  computer  time 
to  generate.  In  retrospect,  it  would  have  been  faster  in 
both  my  time  and  computer  time  to  use  another  second-order 
method  rather  than  continue  with  the  gradient  technique 
for  as  long  as  was  ultimately  required  for  convergence. 

Any  additional  research  on  the  aerodynamic  turn  maneuver 
should  utilize  a  second-order  method  such  as  the  Davidson- 
Fletcher-Powell  variable  metric  technique.  If  continued 
use  of  a  gradient  method  is  desired,  use  of  the  Fletcher- 
Reeves  conjugate  gradient  technique  would  speed  up  convergence 
time . 

Apogee  Analysis 

Before  deciding  to  concentrate  this  study  on  orbits 
with  an  80  km  perigee  a  feasibility  study  of  other  orbits 
was  undertaken.  The  results  of  Harding's  research  indi¬ 
cated  orbits  between  70  and  85  km  perigees  were  profitable 
for  the  aerodynamic  turn  maneuver.  However,  initial  runs 
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at  70  to  75  km  perigee  altitudes  ran  into  difficulty 
when  the  viscous  parameter,  VBAR,  fell  below  the  .01 
value  which  was  the  limit  of  orbiter  aerodynamic 
coefficients  of  lift  and  drag.  Inclination  changes 
at  perigees  of  85  km  and  above  were  insignificant  since 
atmospheric  density  at  these  altitudes  is  too  low  for 
an  effective  aerodynamic  maneuver.  While  it  would  be 
possible  to  compute  lift  and  drag  coefficients  for  perigee 
altitudes  below  70  km,  extended  operations  at  these 
low  altitudes  are  unrealistic  because  of  the  heating 
effects  on  the  orbiter  at  hypersonic  velocities.  This 
tends  to  indicate  the  aerodynamic  turn  maneuver,  while 
feasible,  has  a  limited  range  of  perigee  altitudes  that 
produce  significant  orbit  inclination  changes. 

Comparison  of  the  aerodynamic  turn  maneuver  to  the 
pure  rocket  burn  is  a  purpose  of  this  study.  For  a 
circular  orbit,  the  inclination  change  (Ai)  generated  by 
a  rocket  burn  of  a  specific  change  in  velocity  (Av)  is: 

arcsin  —•  =  Av/2v 

The  value  used  for  v  is  the  scalar  velocity  at  apogee. 
Apogee  is  chosen  since  there  the  shuttle  will  have  the 
lowest  velocity  and  the  largest  Ai  will  be  obtained  for 
a  specified  velocity  change  component  perpendicular  to 
the  plane  of  the  orbit.  Although  the  formula  list  above 
is  for  a  circular  orbit,  the  resulting  Ai  computed  still 


I 


provides  a  good  comparison  value  since  the  shuttle  orbit 
is  nearly  circular. 

The  resulting  Ai  for  a  given  Av  using  a  rocket  burn 
is  listed  in  Table  XI 


Table  XI 


Inclination  Change 


Av  (mps) 

Ai  (deg) 

The  maximum  orbit  plane  change  obtained  by  the  optimized 
controls  in  the  aerodynamic  turn  maneuver  for  the  three 
reference  orbits  with  a  Av  of  60  mps  is  listed  in  Table  XII. 

Table  XII 


Using  Rocket  Motor 


40 

50 

60 

80 

100 

300 

.286 

•  358 

.429 

.572 

.716 

2.15 

Maximum  Orbit  Plane  Change  (Av=60  mps) 


ORBIT 


Ai ( deg) 


80/300  80/500  80/650 


.04 


.44 


.71 


It  is  apparent  that  the  80/300  orbit  accomplishes  almost 
no  inclination  change.  When  compared  to  the  rocket  burn 
results  listed  in  Table  XI  ,  the  80/500  orbit  is  successful 


in  achieving  more  inclination  change  while  the  80/650  orbit 
oetters  the  rocket  burn  by  65%.  The  implications  here  are 
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that  as  the  reference  orbit  gets  more  elliptical  the 
aerodynamic  turn  maneuver  is  more  profitable.  Why  this 
occurs  can  be  explained  by  comparing  the  velocity  histories 
of  the  reference  orbits  inside  the  atmosphere.  As  the 
orbit  becomes  more  elliptical  the  velocity  of  the  shuttle 
increases  in  the  portion  of  the  orbit  near  perigee.  The 
velocity  of  the  shuttle  at  perigee  is  listed  in  Table  XIII 
for  the  reference  orbits. 

Table  XIII 

Velocity  of  Orbiter  at  Perigee 


PERIGEE 
VELOCITY 
(KPS) 

While  this  has  the  effect  of  raising  freestream  Mach 
number,  it  causes  a  reduction  in  the  value  computed  for 
the  viscosity  parameter  (VBAR) .  Table  XIV  lists  the  time 
spent  inside  the  atmosphere  for  each  reference  orbit.  The 
effect  of  the  higher  velocities  in  the  more  elliptical  orbits 
is  quite  apparent  in  the  reduced  time  the  more  ellipitcal 
orbit  allows  the  shuttle  to  stay  in  the  atmosphere. 

Table  XIV 

Time  Inside  Atmosphere 
80/300  80/500  80/650 

TIME 
(sec) 


931 


621 


600 


80/300 


ORBIT 

80/500 


80/650 


7.88 


7-97 


8.00 
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To  understand  why  the  smaller  VBAR  values  result  in 
larger  inclination  changes  it  is  necessary  to  look  at  the 


VBAR  values  calculated  for  the  orbiter.  Table  XV  lists 
the  VBAR  history  for  each  reference  orbit  for  the  portion 
of  the  orbit  between  entry  into  the  atmosphere  and  perigee. 
The  VBAR  history  is  mirrored  for  the  portion  of  the  orbit 
from  perigee  until  it  departs  the  atmosphere. 

Table  XV 

VBAR  History  in  Atmosphere  ( Av=60  mps) 


100 

95 

Altitude 
(km)  :'u 

85 

80 


By  calculating  the  L/D  ratios  for  the  VBAR's  listed 
in  Table  XV  ,  a  comparison  can  be  made  of  the  efficiency 
of  each  orbit.  These  L/D  ratios  are  listed  in  Table  X1/! 


VBAR 


.232 

.098 

.097 

.095 

.094 

.093 

.090 

.086 

.075 

.062 

.062 

.059 

.047 

.045 

OO 

on 

0 

. 

80/300  80/500  80/650 


Table  XVI 

L/D  Ratios  Inside  Atmosphere  (av=60  mps) 
L/D  RATIO  ( a=  25°) 


80/300 

80/500 

80/650 

100 

.761 

1.10 

1.09 

ALTITUDE  9Q 

1.12 

1.13 

1.18 

(km) 

80 

1.35 

1.38 

1.40 

r 


The  L/D  ratios  in  Table  XVI  indicate  that  as  the  orbiter 
enters  the  atmosphere  at  a  higher  velocity  because  of  a 
more  elliptical  orbit,  the  aerodynamic  turn  maneuver  is  more 
efficient.  The  more  efficient  maneuver  obtains  more  inclin¬ 
ation  than  the  less  efficient  maneuver  for  the  same  specified 
velocity  loss. 

At  this  point,  it  is  desirable  to  look  at  the  results 
of  changing  the  acceptable  velocity  loss  per  orbit.  A 
velocity  loss  of  300  mps  is  considered  to  be  the  maximum 
acceptable  Av  loss  in  an  aerodynamic  turn  maneuver.  This 
figure  is  used  since  that  it  the  maximum  Av  the  orbiter 
rocket  motor  can  produce  using  internal  fuel. 

The  inclination  changes  listed  in  Table  XI  which  were 
obtained  using  only  the  rocket  motor  show  that  Ai  is  a 
linear  function  of  Av  for  all  values  of  Av  considered  in  this 
problem.  Based  on  these  results  it  appears  that  accepting 
a  larger  Av  loss  per  orbit  during  the  aerodynamic  turn  maneuver 
would  Rive  a  larger  Ai  proportional  to  the  increase  in  Av. 

Table  XVII  indicates  this  did  not  occur.  In  fact  the 
lower  velocity  loss  returned  a  proportionally  higher  Ai. 

Table  XVII 

Inclination  Change  for  Various  Velocity  Losses  (80/500  orbit) 

VELOCITY  LOSS  (mps) 


30 

40 

50 

60 

Ai(deg) 

.27 

•  35 

.37 

.44 

41 
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These  results  show  that  41  is  not  a  linear  function  of 
4v  for  the  aerodynamic  turn  maneuver. 

When  the  optimization  routine  is  constrained  to 
lose  a  larger  velocity  per  orbit,  there  are  essentially 
two  ways  to  do  this.  One,'  the  orbiter  can  descend  to  a 
lower  altitude  and  stay  longer  in  the  atmosphere.  Since  the 
program  used  here  specified  the  orbit  to  be  optimized, 
this  option  was  not  available.  The  second  method  which 
could  be  used  is  to  increase  the  angle  of  attack.  What 
must  be  analyzed  here  is  the  effect  of  constraints  in 
the  optimization  routine.  Because  the  optimization 
problem  is  set  up  to  obtain  the  maximum  Ai  for  a  given 
Av,  the  program  makes  no  effort  to  find  the  most  efficient 
L/D  ratio  for  each  orbit.  The  optimum  control  history 
for  each  orbit  will  be  applicable  for  a  specific  velocity 
loss.  However,  there  may  be  another  velocity  loss  which 
will  produce  the  most  efficient  Ai  in  terms  of  Av  expended 
to  obtain  the  inclination  change.  Thus,  the  lower  specified 
velocity  loss  returned  a  larger  Ai  on  a  percentage  basis. 

To  understand  why  this  occurs  it  is  necessary  to  calculate 
the  L/D  ratios  at  various  angles  of  attack.  The  L/D 
ratios  for  different  angles  of  attack  at  perigee  VBAR 
are  listed  in  Table  XVIII. 
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Table  XVIII 


L/D  Ratios  versus  Angle  of  Attack  (Perigee  VBAR) 


Angle  of 

Attack 

(deg) 


L/D  RATIOS 


80/300 

80/500 

80/650 

o 

o 

(\J 

1.21 

1.30 

1.39 

30° 

1.19 

1.23 

1.275 

-tr 

O 

O 

•  963 

.982 

1.00 

50° 

.  726 

•  732 

.  794 

Table  X VIE  shows  that  as  angle  of  attack  increases  above 
20°  the  orbiter  L/D  ratio  becomes  less  efficient.  Aero¬ 
dynamic  lift  and  drag  coefficients  are  available  only 
for  angles  of  attack  between  20°  and  50°.  Thus,  the 
most  efficient  aerodynamic  turn  maneuver  for  any  orbit 
in  terms  of  L/D  ratios  would  seek  to  maximize  Ai  for  Av 
which  called  for  an  angle  of  attack  of  20°.  The  optim¬ 
ization  problem  can  be  changed  slightly  so  that  angle 
of  attack  is  constrained  to  remain  at  20°.  Various 
acceptable  Av's  could  be  tried  until  an  optimum  control 
history  for  angle  of  bank  is  determined  that  maximizes 
inclination  change. 

The  inclination  change  which  this  procedure  would 
generate  is  going  to  be  small  (approximately  .20°). 
However,  this  inclination  change  will  be  the  largest  one 
possible  for  the  most  efficient  orbiter  L/D  ratio.  The 
orbiter  can  now  utilize  this  optimum  cont  -ol  history  over 
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several  orbits  until  the  desired  orbit  plane  change  is 

obtained.  Using  this  approach  until  the  total  Av  of 

the  orbits  reaches  the  300  maximum  allowed,  a  large  total 

Ai  is  possible.  This  technique  on  the  80/500  orbit 

produces  a  total  inclination  change  of  3°  while  the  80/650 
.  »  o 

orbit  achieves  a  4  inclination  change.  Thus,  the  Ai 
obtained  by  the  aerodynamic  turn  maneuver  exceeds  the  inclin¬ 
ation  change  of  2.15°  generated  by  the  rocket  burn  (Table  XI) 
by  40%  on  the  80/500  orbit  and  by  86%  on  the  80/650  orbit. 

Higher-Order  Control  Analysis 

The  optimization  program  was  started  with  the  "best 
guess"  of  the  controls  coefficients  as  explained  in  Chapter 
IV  .  The  control  coefficients  were  initially  constant 
values.  The  program  was  allowed  to  optimize  until  the 
increases  in  inclination  change  were  insignificant.  Controls 
were  then  increased  to  three  parameters  for  both  the 
angle  of  attack  and  angle  of  bank  and  the  process  repeated. 

This  procedure  was  accomplished  again  for  five, 
seven,  and  nine  control  parameters  for  only  the  angle 
of  bank.  Higher  order  controls  for  the  angle  of  attack 
did  not  prove  necessary  since  angle  of  attack  essentially 
remained  constant  between  21°-24°  depending  on  the 
reference  orbit. 

The  increases  in  control  parameters  for  the  angle 
of  bank  proved  beneficial  with  each  increase  raising  the 
inclination  change  obtained.  This  occured  through  seven 
control  parameters  for  the  angle  of  bank.  Beyond  seven 


parameters,  there  was  no  appreciable  increase  in  inclination 
change.  The  plots  of  control  history  and  inclination  change 
are  listed  in  Figures  3,4,5,  and  6.  The  increase  in  as  the 
control  parameters  increased  from  one  to  seven  is  approximate! 
.11°.  This  indicates  a  nearly  30%  improvement  in  inclination 
change  can  be  gained  from  the  use  of  higher  order  controls 
rather  than  a  constant  bank  angle. 

The  control  history  for  the  seven  bank  angle  parameters 
indicates  a  converged  solution.  The  oscillations  in  the  bank 
angle  history  for  that  portion  of  the  turn  inside  the  atmosphe 
(.3  <Time  <.6)  have  damped  out.  Further  increases  in  the 
number  of  bank  angle  parameters  had  little  effect  on  the 
control  history. 

An  analysis  of  the  control  history  indicates  the 
converged  solution  commands  the  shuttle  to  roll  inverted, 

<j)  =  l80°,  as  the  atmosphere  is  approached.  This  will  orient 
the  lift  vector  to  aid  in  reaching  perigee.  The  bank  angle 
is  slowly  reduced  until  the  atmosphere  is  reached  with 
4>  =  60  ,  (TIME-.  3).  Most  of  the  inclination  change  occurs 
between  .3<TIME<.6,  where  the  bank  angle  slowly  decreases 
from  80°  to  60°.  As  the  atmosphere  is  departed  (TIME-. 6), 
the  lift  vector  is  oriented  up  ( =  0 0 )  to  assist  in  returning 
to  the  specified  final  altitude.  Flying  this  control 
schedule  is  easily  within  the  capabilities  of  the  orbiter's 


computer . 


This  study  has  shown  that  the  aerodynamic  turn  maneuver 
is  profitable  for  the  Space  Shuttle  Orbiter  and  a  flyable 
control  history  does  result  from  the  use  of  higher  order 
controls  and  optimization.  However,  the  small  inclination 
changes  available  from  the  maneuver  and  the  low  perigee, 
high  eccentricity  orbits  which  must  be  used  to  get  meaningful 
inclination  changes,  severely  limit  the  usefulness  of  this 
maneuver.  Additional  research  on  this  subject  is  not 
recommended. 
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Figure  5*  Control  History-  Five  Control  Parameters 


\ 


H9 


80*0  00 ’Cg 


Bibliography 


1.  Abramowitz,  Milton  and  I.  Segun.  Handbook  of  Mathematical 

Functions .  New  York:  Dover  Publishing,  Inc.,  1964. 

2.  Bate,  Roger,  D.  Mueller  and  J.  White.  Fundamentals  of 

Astrodynamics .  New  York:  Dover  Publishing,  Inc., 
1971. 

3-  Citron,  Stephen.  Elements  of  Optimal  Control.  New  York: 
Holt,  Rinehart  and  Winston,  Inc.,  1969* 

4.  Fox,  Richard  L.  Optimization  Methods  For  Engineering 

Design .  Reading,  Massachusetts:  Addison-Wesley 
Publishing  Company,  1971. 

5.  Harding,  Roger  R. ,  "Optimum  Orbit  Plane  Change  Using 

A  Skip  Reentry  Trajectory  For  The  Space  Shuttle 
Orbiter,"  MS  Thesis,  AF  Institute  of  Technology, 
Dayton,  OH,  December,  1978. 

6.  Hill,  Oliver.  Orbiter  Aerodynamic  Data.  Houston:  1978. 

7.  Hull,  David.  "Application  of  Parameter  Optimization 

Methods  to  Trajectory  Optimization,"  American 
Institute  of  Aeronautics  and  Astronautics,  74-825: 
1974. 

8.  Hull,  David  and  L.  Edgeman.  "Suboptimal  Control  Using 

a  Second-Order  Parameter  Optimization  Method," 
Journal  of  Optimization  Theory  and  Applications , 
17:481-490 :  1975- 

9.  Kuethe,  A.  and  J..  Schetzer.  Foundations  of  Aerodynamics. 

New  York:  John  Wiley  and  Sons,  Inc.,  1959. 

10.  National  Aeronautics  and  Space  Administration.  Space 

Shuttle .  Houston,  1975. 

11.  National  Aeronautics  and  Space  Administration,  United 

States  Air  Force,  United  States  Weather  Bureau. 

U.S.  Standard  Atmosphere,  1962 .  Washington: 

U . S .  Government  Printing  Office,  1962. 

12.  National  Oceanic  and  Atmospheric  Administration,  National 

Aeronautics  and  Space  Administration,  United  States 
Air  Force,  U.S.  Standard  Atmosphere,  1976 .  Washing¬ 
ton:  U.S.  Government  Printing  Office,  1976. 


51 


APPENDIX  A 


Lockheed  Bridging  Formula* 


The  Lockheed  Bridging  Formula  is  used  to  evaluate  the 
aerodynamic  force  coefficients  for  the  Space  Shuttle  Orbiter. 
The  formula  bridges  the  transitional  flow  regime  from 
continuum  flow  to  free-molecular  flow.  The  formula  is: 


trans 


C  .  + 
cont 


(C 


F.M. 


-  C  .  ) 
cont 


•  n ,  ,  . 
sm  (tt  (a 


+  B 


1°s10Kn')  ^ 


where 


C 

cont 

CF.M. 


n 

A 

B 


viscous  force  coefficient  values  at  VBAR  =  0.08 

viscous  force  coefficient  values  at  VBAR  =5.2 

(free-molecular  flow) 


2 

3/8 

1/8 


K  =  Knudsen  number  =  X/L  _ 
n  ref 

Lref  =  12.059  meters 

X  =  mean  free  path  =  RT/P  (2Na^)”^ 

R  =  universal  gas  constant  =  8.312*  x  10  N-m/kg  K 

T  =  temperature 

P  =  pressure 

N  =  Avogadro's  number  =  6.022  x  10^  kmol-'1' 
a  =  effective  molecular  collision  diameter 
=  3-65  x  10-10  meters 
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APPENDIX  B 


Optimization  Program 


5H 


The  program  used  to  solve  the  optimization  of  the 
orbiter  aerodynamic  turn  maneuver  is  listed  in  this  appen¬ 
dix.  A  flow  chart  is  included  to  show  the  calling  order 
of  the  subroutines  (Fig  7).  The  program,  as  listed  here, 
is  set  up  to  solve  the  optimization  problem  using  the 
gradient  technique.  Comments  are  incorporated  into  the 
program  to  explain  the  function  of  the  subroutines  and 
define  key  variables. 


ure  7,  Optimization  Program  Flowchart 
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**♦**»**»  THE  PROGRAM  HZ' IN  COMPUTES  THE  PARTIAL 
DERIVITIVES  USED  IN  THE  PARAMETER  OPTIMISATION  SCHEIE. 

A  SECOND-ORDER  PARAMETER  OPTIMIZATION  TECHNIQUE  IS  JSED 
BASED  ON  THE  PAPER  BY  EDGE MAN  AND  HULL  BY  THE  SAME  NAME* 

******  *  *  * 

PROGRAM  SHUTTLE  (INPUT,  OUTPUT,  TAPES  =INPUT,  PLOT) 

EXTERNAL  F 

DIMENSION  COV(S,  12)  ,CLV(  6,12)  ,ARAY  C2?,5)  ,7B  (6)  »  Yl<6! 
DIMENSION  XI  (6)  ,  YP(6)  ,  WORK (25 J)  ,IWORK( >  ),MNM(2)  ,Y(6I 
DIMENSION  GB<13)  ,GBB(  1G,  1U)  ,  FBT  (10  ) ,  FB3  (1C  ,10)  ,  FBBK10,  10) 

01  ME  NS  I  ON'  03(10)  ,  V  ( lu  )  ,Z  (10) 

DIMENSION  XP(6)  ,  D9P(1C)  ,  RG(4 )  »HBH3T( 4,4 ) j^BHBTI  (4,4i 
DIMENSION  XPP(S) ,XP1<6), XP2( 6) , XP3 C&) , XPA (6) 

DIMENSION  HP (2)  ,HPP(2)  ,HP1  (2  ) ,  HP2  (  2)  ,HP3(?)  ,HP4(2) 

01  MENS  T  ON  H(2),HP(2,lC),HEB(2,lO,10),E(  2,10), 

FI (4,4) , FI ( 4,h) , G  (2) 

DIMENSION  RP(2),P.Q(2)  ,DC  (2)  ,G(1Q),C(2),H8T(10,2) 

DIMENSION  DELINC (7) 

DIMENSION  A  OB  (50)  ,INKCHG  (BO)  ,RADIAN(5Q)  ,AOA(50)  ,  TIME  (50) 
COMMON/B/D(16) 

COMMON /BP/BP (16) 

C0MM0N/N/NE,NP,NPP1,NPP2 

COMMON  U,  RHO,  VISC,  REN,  VPAR,  PHI,  COE  FDP.C-,  COEFLFT 
COMMON  CDV,CLV,ARAY,V1,RA,DT 
COMMON  NOM 
REAL  NOM , INKCHG 
READ*, (XI (I) ,1=1,6) 

INITIAL  CONDITIONS,  XKl  THRU  3)  IS  X,Y,7  POSITION 
INITIAL  CONDITIONS,  XI<4  THRU  6)  IS  VX,VY,VZ  VELOCITY 
RE  AO*,  (  (CDV  (N,M)  ,N=1,  6) ,  M=l,  13) 

RE AD*, ( (CL V (N , M) ,N=1, 8), M=l, 13) 

READ*, ((ARAY(I, J) , J  =  l,5)  ,1=1,22) 

COV  IS  the  DRAG  COEFFiriENT  OATA 
CLV  IS  THE  LIFT  COEFFICIENT  DATA 
ARAY  IS  THE  ATMOSPHERIC  MODEL  DATA 
NE=NUM=»ER  OF  STATE  EQUATIONS 
NP=NUM9ER  OF  PARAMETERS 
NC=NUMBER  OF  CONSTRAINTS 

NPP1  AND  NPP2  ARE  THE  NUM3ER  OF  CONTROL  PARAMETER 
COEFFICIENTS  FOR  ANGLE  OF  ATTACK  AND  BANK, RESPECTIVELY 
N=6 
NE=5 
NP=7 
NC=2 
NPP1=3 
NPP2=3 

C  Cl  AND  C2  ARE  THE  CONSTRAINTS.  HERE  FINAL  ALTITUDE  AND  VEL. 

Cl=6533. 0 
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C 


C2=((XI(t<  ) **2+Xl  (5)**2+XI(6)**2>**  .5)-.  36 
IT=ITI=i 

9  VALUE’S  ARE  COEFFICIENTS  OF  CONTROL  p  ARf  METERS 
8(  i)  =  1  3^1 .456542  23 
B(  2)=23.i<'48603169 
8(  3)=.  00727844741:891 
0(4) =.3104410966147 
B(5)=.3323186 73716 
B(6)=-l. 16885024800 
8( 7) =• 0  0  3^7977  98  997  2 
c  C»S  ARE  INITIAL  LAMBDA  VALUES. 

C<1)=.  OjOi 
C( 2) =• 001 
PC  T  =  0 . 1 

C  OELB  IS  PERTU9ATI0N  STEP  SIZE 

OELB=l . 9  E-C  2 

C  DT  IS  THE  TIME  STEP  SIZE 

DT=.6666565E-03 
ICO  CONTINJE 

PRINT  90, IT 

80  FORMAT (1H,*ITERATI0N  ♦ ,1 3> 

DO  71  1=1, NP 
71  BP  (I )  =  3  ( I ) 

NEXT  STATEMENT  CONTROLS  WHICH  ITERATION  GETS  PLOTTED. 

PLOT  ROUTINE  PLOTS  A  LE  FT  Y  AXIS  (ANGLE  OF  BANK)  VEFSJS 
THE  X  AXIS  (TIME)  AND  A  RIGHT  Y  A  XI 9  (INCLINATION  CHANGE) 
PLOT  IS  6  INCHES  HIGH  BY  5  INCHES  WIDE, 

IF  (IT.  NE.  3)  GO  TO  732 
T  =  0 . 0 

DO  632  N=1 , 23 
INKCHG(N)=ABS(INKCHG(N)) 

632  CONTINUE 

DO  733  N  =  1 , 21 
CALL  THETA(T,UNM) 

A0B(N)=UNM(2)*57 , 29 57 7 95 i 
TIME (N) =T 
T=T*9. 05 
733  CONTINUE 

CALL  PLOT (0  .,-.5 , -3) 

CALL  PLOT  (2.5, 2. 5,-3) 

CALL  PLOT  (5. ,0., -2) 

CALL  PLOT (0 .,7., -2) 

CALL  PLOT (-6., 0., -2) 

CALL  PLOT (0 .,-7. , -2) 

CALL  “LOT ( .5, .5, -3) 

CALL  SSftLE(TXME«R.»21»l) 

CALL  SCME(A0B,6.  ,21,  1) 

CALL  AXIS  (0.,0.,4HTIHE,-4,5. ,3.0, 0.G ,.2) 


I 
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CALL  AXIS(0»»0», 17 MANGLE  OF  BANK  (  B)  , 1? ,6. , 90. , 0 . 0 ,3C « ) 

TIME<22)=0.0 

TIME (23) =.2 

A03(22) =n,C 

A03  (  23)  =  25  , 

CALL  LINE(TIME,AGB, 21,1, 2,29) 

IS=IS-1 

CALL  SCALE(INKCHG,6.,IS, 1) 

CALL  AXIS  <5., 0.,  22HINCL1  NATION  CHANGE  (  »)  , -2  2 ,6  •  ,  91.  P  , 
*INKCHG(IS  +  1) ,INKCHG(IS+2)  ) 

CALL  LINE (T IKE ,1 NKCHG , IS , 1, 2,17) 

CALL  PLOTE(N) 

GO  TO  SP 

732  PRINT  301,3 

801  FORMAT (1X,6E21#12) 

DO  12  1=1, NE 

12  Y1 (I) =X I ( I) 

NOM=O.D 

CALL  EOM(T, Y1,DEG,INKCHG ,IS, -1.0) 

NOM=NOM+i.O 

H(1)=(RA-C1)/10G0. 

H(2)  =  (  <Y1(4)**2+Y1<5)  **2+Yi<6>**2)  ♦*‘.5)  -C2 
PRINT*  ,  "Hill  ANO  H<2)=  ,,,H(1),H(2) 

DELI  NC  ( 1 )  =  -  DEG 
DO  78  K=1,NP 
OO  75  L=1 , NP 
75  9P(L)  =  3(L) 

OBP(K)  =0FLB*3P(K) 

IF<ABS(03P(K)  )  «L  E.DELB)  OBP(K)=DELB 
BP  <K)=B3 (K) +OBP ( K) 

00  13  1=1, NE 

13  XP(X)sXKI) 

SP=2.0 

CALL  EOM(T»XP,DEG,INKCHG,IIP»SP) 

HP (i)=RA-Cl 

HP<2)  =  (<XP<4)**2+XP<5)**2+XPC6)**2>**.5)-C2 
DELINC  C2) =  -D£G 
00  750  L= 1, NP 
753  BP(L)=3(L> 

BP (K)  =  3? ( K) -OBP(K) 

00  14  1  =  1, NE 

14  XPP(I)=XI(I) 

CALL  EOM(T,XPP,DEG, IN KCH G , IIP, SP) 

HPP (1) =RA -C 1 

HPP(2)=<  (XPP(4)**2  +  XPP(5  )  **  2+XPP  CB  )  *  *2)  **.5) -02 
OELINC (3) =-DEG 
00  77  L=1 , NC 

77  :;B(L,K)  =  CHP  (L)-HPP(L)  )/(  2»0*D9P(K)  ) 
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'G8<K)  =  (0ELINC<2>-0ELINC<3)>/  (2.C*D3P(K)  ) 

78  CONTINUE 

C  THIS  SECTION  IS  THE  PARAMETER  OPT IHI7A TIC N  SCHEME 

PRINT*  , **  “ 

PRINT*,”****  H8  MATRIX  ****** 

PRINT*,  (  <H3(I,M>  ,M=i,NP)  ,I=i,NC) 

PRINT*,”  ” 

PRINT*,  •****♦  G8  VECTOR  ****** 

PRINT*,  (  G3<  K)  ,  K=  1 ,  NP) 

PRINT*,"  ” 

00  4  1=1, NP 
00  4  J  =  1 ,  NC 

4  H8T(I,J)=HB(J,I) 

DO  5  I  =  1 ,  NC 

DO  5  K  =  1 , NC 
H3H8T ( I , K )  =  C  •  0 
DO  5  J=l, NP 

5  H8H3T ( I ,  K ) =  HBHBT <I,K) +HB  (1, J ) *HBT ( J, K) 

CALL  GAUSO(NC,i« C E- 30 , HG HBT , HBHBTI ,KER, 4) 

00  655  1=1, NC 

RG(I)=0. C 
DO  665  J=1 , NP 

666  kG(I)  =  RG(I)+HB(I,J)*G8(J) 

SS=.90OOOCi 
0=.001 
00  7  I  =1 , NC 
C<I)=0.0 
00  7  J=1,NC 

7  C(I)=C(I)  +HBH3TI  (I,J)*<0*H(J)  -RG  ( J  ) ) 

11  CONTINUE 

PRINT*,"  " 

PRINT*,  ”***♦  C  VECTOR  ****** 

PRINT*,"  " 

PRINT*,"  **,C<i>,"  ”»C  (2) 

00  19  1  =  1, NP 
CC=0.0 

DO  110  J=i, NC 
110  CC=CC*C(J)*H9(J,I) 

10  F3T (I ) =GP (I ) +CC 
PRINT*,"  " 

PRINT* ,”*********** F9T  VECTOR  ******  ******  ********* 
PRINT*,"  " 

PRINT*, (FBT (K) ,«=1,NP) 

00  549  1=1, NP 
00  549  J=l, NP 
F88I (I , J) =0  •  9 
IF(I.EQ.J)  FB8I (I,J)=1«0 
549  CO NT IN JE 
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00  20  1  =  1, NC 
00  20  <  =  1, NP 
E(I»  K)  =0.'i 
00  20  J  =  i  ,NP 

20  E(I,<)  =  E<I,K)  fFBRI  (J,K)*HE(I,J) 

00  2  5  1=1, NC 
00  25  K=i,NC 
Fi(i,o=r.o 
00  25  J=1,NP 

25  Fl(I»K>=Fi(I*K)  +  E<I,J)*HBT(J,K) 

CALL  G AUSD ( NC , 1*  0  E-30 , Fi , FI ,KER»  4) 
00  33  1=1, NC 
GCI) =3 • C 
00  30  J  =  1 , N P 

30  GCI)=G(I> +E(I,J) *F3T( J) 

00  546  1=1, NC 

548  DC (I) =0. 0 

00  +0  1=1, NP 
S(I>=0.1 
00  40  J= 1 , NC 

40  S(I)=S(I) +HBT (I, J)*DC (J) 

00  45  1  =  1, NP 
V(I)=C.O 
T(  I) =0  •  0 
00  45  J=1 , N P 

tf(I)=V(I) +FBBI (I,J)*FBT(J) 

TCI)  =?  (I)  +F9BICI,  J)*S(J) 

45  D9(I)=-<SS*V(I))~ZCI> 

O0N=O.O 
00  145  1=1, NP 
145  08N=DBN*09(I)*0B(I) 

03N=SQRT  C09N) 

BN=0 • 0 

00  255  1=1, NP 
8N=9N+9(I) *8 ( I > 

255  CONTINJE 

BN=SORT (BN) 

HN  =  G  »  9 

00  150  1=1, NC 
150  HN=HN+H(I)*H(I) 

HN=SQRT(HN) 

PRINT*,"  " 

PRINT*,  "****  08N  HN  ♦***” 

PRINT*,"  ** 

PRINT  3  2 , 03N,HN 


82  FORMAT (4X, 2E26.8 ,/) 
PRINT* " 
PRINT*.-****  03  ****- 


r 
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PRINT*,**  ** 

PRINT  33,93 

83  FORMAT  ('+X,6E20. 8,/) 

PRINT*,**  ** 

PRINT*,******  OC  ****** 

PRINT*,  **  ** 

PRINT  64,00 

84  FORMAT  UX,  3E23. 3,///) 

IFU03N.LT. 1.3E-04)  .AND.  ( HN. LT. i. 0 E-03)  )  00  TO  60 
DO  50  1=1, NP 
50  9(1) =3(1) +08(1) 

IF  (9(2).  GT  .50.0)  9(2)=50.0 
IF  (9(2). LT*  20.0)  B(2?=2C.O 
IF  (8(3). LT. -3  .14)  B(3)=-3.i4 
IF(9(3)  .GT.  3.14)  B(3)=3.l4 
DO  55  1=1, NC 
55  C(I)=C(I) +DC(I> 

IT  =  IT  +  1 

IF  (IT.  GT  •  ITI  +  25)  GO  TO  60 
GO  TO  100 
60  CONTINUE 

PRINT  83 

85  FORMAT ( ?  X*  CONVERGENCE  *> 

PRINT  35,3 

86  FORMAT (//,4X,6E21«12, //) 

STOP 

ENO 

C  *****SU3ROUTINE  GAUSD  IS  A  MATRIX  INVERSION  ROUTINE***** 

SUBROUTINE  GAUSD (M, EPS, B 9 ,CC ,K£R, LAY ) 

DIMENSION  8B(LAY,LAY) ,CC (LAY ,LAY) , A( 20 , 20) ,X(20,20) 

EP=EPS 

N=M 

DO  130  J=1,N 
OO  100  <  =  1,N 
100  A(J,K)=3B(J,K) 

DO  1  1  =  1,  N 
DO  1  J=1 , N 

1  X(I,J)=0.0 
00  2  K=l, N 

2  X(K,K) =1 ,G 
00  34  L=i,N 
KP=J 

2=3.0 

00  12  K=L,N 

IF(Z-ABS(A(K,L)))11,12,12 
11  Z=  ABS ( A  ( X, L  ) ) 

KP= K 


12  .CONTINUE 


13 


IF(L-KP)  13,20,  20 
00  14  J=L,N 
7=A(L, J) 

A(L, J) =  A<  KP,J) 

14  A(KP,J)=^ 

DO  15  J  =  1 ,  N 
E=X(L, J) 

X(  L  ,  J)  =X(KP,J) 

15  X(KP,J)=7 

2C  IF(A8S(A  (L,L))-EP)5C,£,0,  3G 

30  IF(L-N)31,3i,3<* 

31  LP1=L+1 

DO  36  K=LP1,N 
IF(A(K,L))32»35»32 

32  RATIO=A(K,L)/A(L,L) 

00  33  J=LP1,N 

33  A(K,  J)  =MK,  J)-RATIO«A  (L,  J) 

DO  35  J  =  1 , N 

35  XtK, J)  =X (K, J) -RATIO*X  (L,  J) 

36  CONTINUE 

34  CONTINJF. 

00  43  1=1, N 
II=N+1-I 
DO  43  J= 1  ,N 
S=0  •  9 

IF (II- N) 41,43,43 

41  IIP1=II+1 

00  42  K= IIP1, N 

42  S=S+A(II,K)*X(K, J) 

43  X(II,J)=(X(II,J) -S)/A (II, II) 

KER=1 

00  200  J  =  1,N 
00  200  <=1,N 
200  CC <J,K)=X (J,K) 

GO  TO  75 
50  KER=2 

70  PRINT  71 

71  FORMAT ( IX ,• MATRIX  SINGULAR  IN  GAUSO*) 

75  CONTINUE 

RETURN 

END 

C  ♦•♦♦SUBROUTINE  THETA  COMPUTES  THE  ANGLE  OF  ATTACK  A 'IP  3ANK 

C  ANGLES  USING  THE  CHE3YS  CHEV  POLYNOMIAL  FORM  FROM  SJB.OHEBY. 

SUBROUTINE  THETA (T,UNM) 

01  MENS  ION  UNM ( 2) 

COMMON/ 3 P/BP (16) 

COMMON/N/NE ,NP,NPP1,NPP2 
COMMON/TSH/TOT (8) 
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CALL  CHEBY(T) 

UNM(i)=D.O 
00  20  I=i,NPPi 

20  UNM(i)=JNM(i)+BP(I  +  l)*TOT  (I) 

UNH(2)=0.0 
00  23  1  =  1,  NPP2 

23  UNM ( 2) = JNM  < 2) +9P (I+NPP1+ 1)*T0T(I> 

RETURN 

END 

C  ***** SUBROUTINE  CHEBY  COMPUTES  THE  CHEBYSHEV  P0LYN0-1IALS***** 

SUBROUTINE  CHEBY (T) 

COMMON/TSH/TOT (8) 

TOT ( i) =1  • 0 
TOT(2)=2.f-*T-i.O 
T0T(3)=3.l*T**2-8.C*T+i.  C 
TOT (4) =32.C*T**3-48 .3*T* *  2  + 18. 0*T-i . 1 
TOT  (5)  *126.tf*T**A-256.I)*  T**  3+150 . 3  *T**  2  -32.  U*T+1.0 
T0T(6)=E12.C*T**  5-123C.e*T«Mt4  +  112G.C*T**3-4J3.0*T**2*5J.  O'T-l.O 
TOT  (Z) =2048 •0*T**6-bi44*  C^T+^E  +  BBi  2.  0 *  T  **-3584 . 0*T»  '3 
1 +640. 3 *T** 2-72. 0*T  +  i. 0 
T0T<8) =n.C 
RETURN 
END 

***►♦ SUBROUTINE  EOM  COHFUTES  THE  INCLINATION  CHANGE(DEG). 

IT  USES  THE  CC6E00  ROUTINE  ODE  TO  DIFFERENTIATE  THE 
EQUATIONS  OF  MOTION  WHICH  ARE  CONTAINED  IN  SUBROUTINE  F. 
SUBROUTINE  EOM (T , Y , OE G,I NKCHG, IIP, SP) 

EXTERNAL  F 
COMMON/B/9M.O) 

C0MM0N/BP/3P(16) 

C0MM0N/N/NE,NP,NPPi,NPP2 

COMMON  J, RHO, NISC,REN , VPAR, PHI,CDEFDRG,  COFFLFT 
COMMON  Cnv,CLV,ARAY,Vl,RA,DT 
COMMON  NOW 


2 

3 


DIMENSION  I NKCHG (50 ) 

DIMENSION  CDV(8, 12) ,CLV( 8,12) , ARAY (22,  S), ZB (8) 
OIMENSION  Y(6),YP<6),  WORK  (2BO)»IWORK(5)  ,UNM(2) 
REAL  I NKCHG 
IIP=1 
TTPsO.O 
00  2  1*1,6 
Y(I)*Y(I)*1.0E+G3 
DO  3  1*1,6 
28  (I )  =  Y  ( I ) 

RELERRsl.0E-C6 
T=  0  . 

TOUT=J • 

IFLAG*! 
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KOUNT  =  9 

1  CONTINUE 

TOUT=TOT 

CALL  0DE(F,6,Y,T,T0UT,RELERR,REL£RR,  IPL  AG » WORK,  IHORO 
IF(V3AR.LT.  0.005  )  GO  TO  50 
IF ( NOM-1 • 0  > 48,44 , *4 
46  CONTINUE 

40  Ai=ZB(2>47B(6)-ZB(3)*ZB(£) 

B1=Y (2>  4 Y (6 )  -Y (3 )  4Y  (5) 

Cl=Z8( 3)  479 (M -ZB(1)*ZB(6) 

01=Y(3V4Y(4)-Y(i)4Y(6) 

El  =  Z9(i)478(5>-ZB<2>4Z8(4) 

Fi=Y(i)fY(5)-Y(2)4Y  (4) 

G=Al*'3ifCi*Di*El*  FI 
H=(Ai4»24Ci442  +  Ei442>  440  .5 
0=(Bi4424Di442  +  Fi442) 44C  .5 
CI=G/( H40) 

IF(AB3(CI>.GT.i.  )  CI=i. 

AI=ACO  3  ( C  I) 

OEG=AI4180. 0/3.14159 
IF  (SP.  GT  .0.0)  GO  TO  736 
IF (TOUT. GE. TTP)  GO  TO  735 
GO  TO  736 

735  CONTINUE 
INKCHG  (IIP)  =OEG 
TTP=TTP«-.G5 
IIP=II3*1 

736  DO  4  <=1,6 

4  Y(K)=Y(K)41.0E-03 

VEL=(Y(4)*4  2+Y (5)V424Y(6)4*2)44.5 
IF (TOUT. GE. 1. 0 )  GO  TO  49 
IF (KOUNT • GT  •  0 )  GO  TO  45 
PRINT4, “THIS  IS  NOMINAL  DATA” 

GO  TO  46 

45  IF (KOUNT .LT.50)GO  TO  44 

46  CONTINUE 
APER=U/1000. 

IF(KOUNT.GE.5Q>  KOUNT=i 

44  KOUNT= KOUNT +1 

IF (TOUT . L E . 1. 0)  GO  TO  1 
IF (TOUT. GE. 1, 0 )  GO  TO  4f 

49  PRINT4, “THE  GEOCENTRIC  ALT  IS  ",RA 
RETURN 

50  PRINT4, “V^AR  LESS  THAN  C.liGS  ”,VBAR 
RETURN 

END 

C  *♦»*#  SUBROUTINE  F  CONTAINS  THE  EQUATIONS  OF  MOTION  WHICH 

C  ARE  USED  BY  ODE. 


SUBROUTINE  F(T,Y,YP) 

COMMON  U,PHO,VISC,REN,VeAR,PHI,CO£FORG,  COFFLFT 
COMMON  :OV,CLV,ARAY,Vi,F A,OT 
COMMONA3P/8P(1S) 

DIMENSION  Y  (6)  ,YP(6),CDV  (8,12)  ,CLV  (E,  ,12  )  ,  ARA  YC22 , 5) ,  UNM  (  2 ) 

GM=3.98f 2216E+14 

OM=7.292ii5BE-05 

SA=9d. 3 

SM=975  0 " «  C 

RA=(Y(1)**2+Y(2>«*2+Y<3>  **2>**.5 
VS  =  YU)  **2  +  Y<5>**  2+YIC,)'  *2 
V1=VS**,5 

CALL  ATMOS (TM, WMOL) 

CALL  VISC°AR(TM, WMOL) 

CALL  THET A ( T ,UNM ) 

CALL  BIVINT (ALPHA, UNM) 

PHI=UNM(2) 

SPHI=3IN(°HI) 

CPHl=COS(PHI) 

SLIFT= • 5*COEFLFT*  SA*RHO*  VS 
DRAG=» 5* COE  FORGES A* RHO* VS 
SLM=SL  I FT  /SM 
SOM=ORA3/SM 
AV=Y(4)-OM*Y(2) 

BV=Y (5  )  +0  M*  Y  ( 1 ) 

CV=Y(S) 

AH=BV*Y(3)-CV*Y( 2) 

BH=-(AV*Y(3)-CV*Y(1)) 

CH=AV* Y(2)-BV*Y(i) 

AL=CV*3H-BV*CH 
BL=-(CV*  AH-CH*  AV) 

CL=BV*AM-BH*AV 
SV=(AV*»2+3V**2+CV**2)**  .5 
SH=( AH**  2+8H** 2+ CH**2) ** .5 
SL=(AL**2  +  BL**2+CL**2)**  .5 
R3=(Y(1)**2+Y(2)**2+Y (3)**2>**1.5 
YP(1)=Y(A)*BP<1) 

YP(2)*Y(5)*BP(1> 

YP(3)*Y(6)*BP(i) 

YS4=-(GM/R3*Y(1) ) +SLM* <-SFHI*AH/SH+CPHI*AL/SL)-SDM*\V/SV 
YP(V)sYS4*3P(l) 

YS5*- ( GM/P3*Y ( 2) ) +SLM* (-SPHI*3H/5H fCOHI *RL/SL> -SDM*3 V/3 V 
YP(5) *Y35*8P<1) 

TS6=-(GM/R3*Y(3) ) +SLM* (- SPHI*CH/SH  +C PHI "CL/SL) -SDM*: V/3 V 
YP(S)sYS5“BP(l) 

RETURN 

ENO 
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irv  a)  <0 


r 


c  ♦»*** SUBROUTINE  ATMOS  COMPUTES  THE  MOL  EC. WEIGHT ( WMOL )  A. 

C  OENISTY  OF  THE  ATMOSPHERE 

SU3R0'JTINE  ATMOS  < TH ,WMOL  ) 

COMMON  U,PHO, VISC, REN, VBAR, PHI, COEFDRG,  COEFLFT 

COMMON  CDV,CLV,ARAY,V1»RA,DT 

DIMENSION  COV(  8,12)  ,CLV(  6  ,12)  ,ARAY  <22,3  ) 

GM=3.9BS2216E+14 
R=8.31432f+C3 
RH0Z=1 .2250 
TM7*  238.15 
U=RA-635  6766. 0 
00  7  M=1 , 22 
IF (U-ARAY (M,l) )  5 ,6,7 
7  CONTINUE 

M=23 
M=M-1 

IF<ARAY<M,4)) 8,9,8 

TM=ARAY(M,2)+ARAY(M,4>  *<U-ARAY<M,1)) 

QsGM*ARAY  (M,5)  /(R*RA**3) 

SI GMA=  EXP ( ( 1.0  +  (C/ARAY(M,4) ) ) * < ALOG < ARA Y(M, 2 )/TM) )  )*£RAY 
-3) 

GO  TO  10 

9  TM=ARA Y ( M, 2 ) 

0=GM*ARAY<M,5)  /<R*RA*  +  2) 

SIGMA=ARAY(M,3)*£XP(-<Q*  (U-ARAY  (M,  1)  ) )  /  ARA Y  ( M, 2)  ) 

10  RH0=RH07*siGMA 
NMOL=ARAY (M ,5 ) 

RETURN 

ENO 

C  *****SJRROUTINE  VISCPAR  COMPUTES  VB AR***B* 

SUBROUTINE  VISCPAMTM, WMOL) 

COMMON  U,RHO, VISC, REN, VEAR, PHI, COEFORG, COEFLFT 

COMMON  3DV,CLV,AF.AY,Vi»KA»DT 

DIMENSION  CDV (8,12) ,CLV< 8,12) ,ARAY (22,5) 

T=TM*WM0L/28.9644 

VISC=(1.458E-06*T**1.&)/  (T+110.4) 

REN*RH  0*  VI*  32. 77/ VISC 

A=SQRT (1. 15*8»31432E+L 3* TM/ 28.9644) 

TP*.458*T  +  1366.3+0.01462&*<  < Vl/A)**2) 

EXPl=-5.0/T 

EXP2=-5.C/TP 

ANUM=T  +122.1*10**  EX  PI 

AOEN=TP*122.i*15**EXP2 

CP*((TP/T)**.5)*ANUM/A0EN 

V9AR=( Vl/A)  *(  CP/REN)**  .5 

RETURN 

END 


C  ROUTINE  BIVINT  IS  THE  BI VARIANT  INTERPOLATf ON  SCHEM 

C  USED  TD  GENERATED  THE  ORPITER  COEFFICIENTS  OF  LIFT  ANO  DRAG 

SUBROUTINE  BIVINT (ALPHA, UNM) 

COMMON  U,PHO, VISC, REN  jVEAR, PHI, COEFORG , COcFLFT 
COMMON  OOV»CLV»ARAY,V1»F'A,QT 

DIMENSION  COV(8,12)  ,CLV(  8,12)  ,ARAY  (2  2, 5  )  ,UNM  (2) 

IFCVBAR.GT.5.2)  GO  TO  1C 

IF (UNM (1) .LT.20. C)  UNH(1)=2Q.C 

IF  (UNM (1)  .GT.50.C)  UNM  (1  )  =50  .0 

ALPHA=  UNM ( 1 ) 

00  3  1=2,8 

IF  (ALPMA-CDVd  ,1)  )  2,2 ,3 

3  CONTINUE 

2  DA =C0V(I,1) -ALPHA 

DO  5  J=? , 12 

IF (VBAR-COV (1, J) )4,4, 5 

5  CONTINUE 

4  DV=CDV(1,J)-V3AR 
OELV=COV (1, J)-COV(i, J-l> 

P=OV/OEL  V 

Q=0A/5 • 3 

COEFDRGs  ( 1-P) *  (1-Q)  *CDV(  I  ,J)  +PM1-Q)  *CD  V(I,  J-l>  +QM1-P)  * 
-CDV(I-1,  J>+P*Q*CDV(I-i,J-l> 

00  7  1  =  2,8 

IF(ALPMA-CLV(I,1))6,6,7 

7  CONTINUE 

6  DA=CLV ( T , 1 ) -ALPHA 
DO  9  J=2 , 12 

IF ( VBAR-CLV (1, J) ) 8,8,9 

9  CONTINUE 

8  OV=CLV(l, J) -VBAR 
0ELV=CLV(1,  J)-CLV(1,J-1) 

P=DV/OELV 

Q=  OA/5 . 3 

COEFLFT  =(i-P)*(i-Q)*CLV(I,J)  +P*  (i-QJ’CLVd,  J-l)  +Q*  (l-P)  *CLV 
-(I-1,J)4-P*Q»CLV<I-1,J-1) 

GO  TO  11 

10  COEFDR3=0.0 
COEFLFT  =  0  «  0 

11  RETURN 
ENO 
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